// Diffractor.cpp: implementation of the Diffractor class.
//
//////////////////////////////////////////////////////////////////////

#define PI   3.1415926535897932385

#include "Diffractor.h"

#include <fstream>  // used to get ifstream
#include <sstream>
#include <algorithm>

//extern "C"
//{
#include "../Atominfo/atominfo.h"
//}

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////


Diffractor::Diffractor(char *DiffractorFile,
					   Position *W1, Position *W2, Position *W3, char **atomLabel, int nAtom)
{

  NAtom = nAtom;  // note NAtom is redefined in Diffractor_read_old_style()

  string filetype(DiffractorFile);

  size_t found = filetype.find(".hcv_dash");
  if (found != string::npos)
    Diffractor::Diffractor_read_hcv(DiffractorFile, W1, W2, W3, atomLabel);
  else
    Diffractor::Diffractor_read_old_style(DiffractorFile, W1, W2, W3, atomLabel);
}




Diffractor::~Diffractor()
{
delete  Data;

deleteA(Intensity);
deleteA(ReAmplitude);
deleteA(ImAmplitude);
deleteA(H);
deleteA(K);
deleteA(L);
deleteA(DH);
deleteA(DK);
deleteA(DL);
deleteM(FiniteFactor);
deleteM(FDeriv);

deleteSymmA(SinX, MaxH);
deleteSymmA(CosX, MaxH);
deleteSymmA(SinY, MaxK);
deleteSymmA(CosY, MaxK);
deleteSymmA(SinZ, MaxL);
deleteSymmA(CosZ, MaxL);
}


//////////////////////////////////////////////////////////////////////
// functions
//////////////////////////////////////////////////////////////////////

double *Diffractor::Build(Position *AtomPosition)
{
FourierTransform(AtomPosition);

for(int n=1;n<=NIntensity;n++)
	Intensity[n] = ReAmplitude[n]*ReAmplitude[n] + ImAmplitude[n]*ImAmplitude[n];

return (Intensity);
}


////////////////////////////////////////
// private methods
////////////////////////////////////////

void Diffractor::Diffractor_read_hcv(char *DiffractorFile,
					   Position *W1, Position *W2, Position *W3, char **atomLabel)
{

  // See comment top for why this code commented out

  cout << "\nDiffractor_read_hcv method code commented out for now\n";
  exit(1);




  // Number of off-diagonal weight matrix elements read in for each line
  // of the 'new' style .hcv file.
  // Notice when a multiplet SPLIT it into singlets
  vector<int> num_cor_element;




  // read diffraction data into attribute data which has
  // the structure:

  int*    h;                   // pointer to an array of h's
  int*    k;                     
  int*    l;
  double* ave_intensity;       // pointer to an array of the likelihood expectation values of the intensities
  double **inv_cov;            // pointer to inverse of the likelihood covariance matrix of the intensities
  int N;   // the number of intensities





  populate_store_data(DiffractorFile, num_cor_element);


  FracX = (2.0*PI)*(*W1);
  FracY = (2.0*PI)*(*W2);
  FracZ = (2.0*PI)*(*W3);

  M = TPMatrix(FracX,FracY,FracZ);

  NIntensity = data.N;

  Data = new DataHolder(NIntensity);

  newA(Intensity,		NIntensity);
  newA(H,				NIntensity);
  newA(K,				NIntensity);
  newA(L,				NIntensity);
  newA(DH,			NIntensity);
  newA(DK,			NIntensity);
  newA(DL,			NIntensity);
  newA(ReAmplitude,	NIntensity);
  newA(ImAmplitude,	NIntensity);
  newM(FiniteFactor,	NAtom, NIntensity);
  newM(FDeriv,		NAtom, 3);

  // to calculate average form factor decay

  vector<double> averageDecay;


  double weightMatrixOffDiagonalCutoff = 20.0;
  
  // populate Data.Hessian from inv_cov array
  // as read in from hcv file, where diagonal elements
  // are i/sigma and off diagonal elements percentages

  for (int i = 1; i <= NIntensity; i++)
  {
    Data->Hessian[i][i] = data.inv_cov[i][1]*data.inv_cov[i][1];
  }

  // filling in the off diagonal elements

  for (int i = 1; i <= NIntensity; i++)
  {
    for (int j = i+1; j-i+1 <= num_cor_element[i]+1; j++)
    {  
      if (data.inv_cov[i][j-i+1] > weightMatrixOffDiagonalCutoff)
      {
          // convert percentage inv_cov into absolute hessian
          Data->Hessian[i][j] = data.inv_cov[i][j-i+1]
            * sqrt(Data->Hessian[i][i] * Data->Hessian[j][j]) / 100.0;

          Data->Hessian[j][i] = Data->Hessian[i][j];
      }
    }
  }



  MaxH = 0;
  MaxK = 0;
  MaxL = 0;

  for(int m=1;m<=NIntensity;m++)
	{
	  Data->Average[m] = data.ave_intensity[m];
  	
	  H[m] = data.h[m];
	  K[m] = data.k[m];
	  L[m] = data.l[m]; 


	  DH[m] = static_cast<double>(H[m]);
	  DK[m] = static_cast<double>(K[m]);
	  DL[m] = static_cast<double>(L[m]);



	  if(abs(H[m])>MaxH) MaxH = abs(H[m]);
	  if(abs(K[m])>MaxK) MaxK = abs(K[m]);
	  if(abs(L[m])>MaxL) MaxL = abs(L[m]);
	}

  newSymmA(SinX, MaxH);
  newSymmA(CosX, MaxH);
  newSymmA(SinY, MaxK);
  newSymmA(CosY, MaxK);
  newSymmA(SinZ, MaxL);
  newSymmA(CosZ, MaxL);


  SinX[0] = 0.0;		CosX[0] = 1.0;
  SinY[0] = 0.0;		CosY[0] = 1.0;
  SinZ[0] = 0.0;		CosZ[0] = 1.0;


  // creating your own FiniteFactors

  for(int m=1;m<=NIntensity;m++)
  {
     Position temp_pos = (*W1)*H[m]+(*W2)*K[m]+(*W3)*L[m];
     double stol2 = temp_pos*temp_pos / 4.0;

	   for(int n=1;n<=NAtom;n++)
     {  
        const T_SF_WK95_CAA *AtomScatFac;
        AtomScatFac = FindSF_WK95_CAA(atomLabel[n], 0);
        FiniteFactor[n][m] = CalcSF_WK95_CAA(AtomScatFac, stol2);
     }

     // optionally include a debye-waller fall off if this has not 
     // already been done in hcv file. 
     // In general depending on data this may have little effect
     // The number 2.42 seems a bit arbitray and of course will
     // be heavily influenced by temperature
	   for(int n=1;n<=NAtom;n++)
     {  
       //FiniteFactor[n][m] *= exp(-2.42*stol2);
     }
  }

  // Make form factors constants by applying mild approximation
  convertFormFactorDecay(atomLabel);


}


void Diffractor::Diffractor_read_old_style(char *DiffractorFile,
					   Position *W1, Position *W2, Position *W3, char **atomLabel)
{
  int    i, j, n, m, NElement, NSpecial;
  double d1, d2, d3, dummy;

  SAFEifstream infile(DiffractorFile);


  FracX = (2.0*PI)*(*W1);
  FracY = (2.0*PI)*(*W2);
  FracZ = (2.0*PI)*(*W3);

  M = TPMatrix(FracX,FracY,FracZ);

  infile >> NIntensity >> NElement >> NAtom >> NSpecial;

  Data = new DataHolder(NIntensity);

  newA(Intensity,		NIntensity);
  newA(H,				NIntensity);
  newA(K,				NIntensity);
  newA(L,				NIntensity);
  newA(DH,			NIntensity);
  newA(DK,			NIntensity);
  newA(DL,			NIntensity);
  newA(ReAmplitude,	NIntensity);
  newA(ImAmplitude,	NIntensity);
  newM(FiniteFactor,	NAtom, NIntensity);
  newM(FDeriv,		NAtom, 3);


  for(m=1;m<=NElement;m++)
	  {
	  infile >> i >> j >> d1 >> d2 >> d3;
  	
	  if(i>NIntensity || j>NIntensity)
		  {
		  cerr << "Problem with Hessian elements in *.inf!\n\n";
		  exit(1);
		  }

	  if(i==j)	
		  Data->Hessian[i][i] = d1;
	  else
		  {
		  Data->Hessian[i][j] = d1/2.0;
		  Data->Hessian[j][i] = d1/2.0;
		  }
	  }



  MaxH = 0;
  MaxK = 0;
  MaxL = 0;

  for(m=1;m<=NIntensity;m++)
	  {
	  infile >> d1 >> d2 >> d3 >> Data->Average[m];
  	
	  //for reading in special .inf
  /*	int int_d1, int_d2, int_d3;
	  infile >> i >> j >> int_d1 >> int_d2 >> int_d3;


	  H[m] = int_d1;
	  K[m] = int_d2;
	  L[m] = int_d3;

	  DH[m] = double(H[m]);
	  DK[m] = double(K[m]);
	  DL[m] = double(L[m]);*/

	  DH[m] = ceil(d1/(2.0*PI) - 0.5);
	  DK[m] = ceil(d2/(2.0*PI) - 0.5);
	  DL[m] = ceil(d3/(2.0*PI) - 0.5);

	  H[m] = int(DH[m]);
	  K[m] = int(DK[m]);
	  L[m] = int(DL[m]); 

	  if(abs(H[m])>MaxH) MaxH = abs(H[m]);
	  if(abs(K[m])>MaxK) MaxK = abs(K[m]);
	  if(abs(L[m])>MaxL) MaxL = abs(L[m]);
	  }


  newSymmA(SinX, MaxH);
  newSymmA(CosX, MaxH);
  newSymmA(SinY, MaxK);
  newSymmA(CosY, MaxK);
  newSymmA(SinZ, MaxL);
  newSymmA(CosZ, MaxL);


  SinX[0] = 0.0;		CosX[0] = 1.0;
  SinY[0] = 0.0;		CosY[0] = 1.0;
  SinZ[0] = 0.0;		CosZ[0] = 1.0;


  for(m=1;m<=NIntensity;m++)
	  for(n=1;n<=NSpecial;n++)
		  infile >> dummy;

  // for reading in special .inf

  /*for(m=1;m<=NIntensity;m++)
	  for(n=1;n<=NSpecial;n++)
		  infile >> fisse;*/

  for(m=1;m<=NIntensity;m++)
	  for(n=1;n<=NAtom;n++)
		  infile >> FiniteFactor[n][m];

  // Make form factors constants by applying mild approximation
  convertFormFactorDecay(atomLabel);
}


/*
Approximately convert form factor decay from the form factors to the
data. The form factor decay will not be exactly the same for every
type of atom - but fairly similar - in particular from the point of
view of structure solution. Take the average form factor decay and
apply to data and transfer the form factor constants as function of
d-spacing. 
*/
void Diffractor::convertFormFactorDecay(char **atomLabel)
{
  /// print out FiniteFactor for curiosity before convertion
/*  ofstream fisse("before_compare_finiteFactor.txt");
	for(int n=1;n<=NAtom;n++)
  {
    for(int m=1;m<=NIntensity;m++)
    {
		  FiniteFactor[n][m] = FiniteFactor[n][1];
		  fisse << FiniteFactor[n][m] << " ";
    }
    fisse << endl;
  }
  fisse.close();*/


  // get form factors at q=0.0
  vector<double> qZero;
  qZero.push_back(0.0); // because all other vectors in this code starts with 1
	for(int n=1;n<=NAtom;n++)
  {  
    const T_SF_WK95_CAA *AtomScatFac;
    AtomScatFac = FindSF_WK95_CAA(atomLabel[n], 0);
    qZero.push_back(CalcSF_WK95_CAA(AtomScatFac, 0.0));
  }

  // rescale form factor 'matrix' so all form factors scale to 1 electron
	for(int n=1;n<=NAtom;n++)
  {
    for(int m=1;m<=NIntensity;m++)
    {
		  FiniteFactor[n][m] /= qZero[n];
    }
  }

  // calculate average form factor decay and store in 1st row
  // of FiniteFactor (i.e. atom 1)
  for(int m=1;m<=NIntensity;m++)
  {
    double average = 0.0;
	  for(int n=1;n<=NAtom;n++)
    {
      average += FiniteFactor[n][m];
    }
    FiniteFactor[1][m] = average / NAtom;
  }

  // In addition I would prefer the form factor constant to be
  // such that the highest form factor is 1.0 since this will
  // ensure that F[n,H] approximate a constant function as n
  // is increased
  double qZeroMax = *max_element(qZero.begin(), qZero.end());
  for(int m=1;m<=NIntensity;m++)
  {
    FiniteFactor[1][m] *= qZeroMax;
  }

  // now transfer form factor decay to data i.e. take
  // W_hk * (I^o_h - a_h*I_h) * (I^o_k - a_k*I_k) and tranform into
  // W_hk*a_h*a_k * (I^o_h/a_h - I_h) * (I^o_k/a_k - I_k)
  for (int i = 1; i <= NIntensity; i++)
  {
    Data->Hessian[i].multiplyElementByElement(FiniteFactor[1]);
    // multiply each row element by FiniteFactor[1][i]
    Data->Hessian[i].multiplyElementByElement(FiniteFactor[1][i]);

    Data->Average[i] /= FiniteFactor[1][i];
  }

  // finally leave form factors 'matrix' as constants
	for(int n=1;n<=NAtom;n++)
  {
    double scale = qZero[n] / qZeroMax;
    for(int m=1;m<=NIntensity;m++)
    {
      FiniteFactor[n][m] = scale;
    }  
  }  


  /// print out FiniteFactor for curiosity before convertion
 /* ofstream fisse2("after_compare_finiteFactor.txt");
	for(int n=1;n<=NAtom;n++)
  {
    for(int m=1;m<=NIntensity;m++)
    {
		  FiniteFactor[n][m] = FiniteFactor[n][1];
		  fisse2 << FiniteFactor[n][m] << " ";
    }
    fisse2 << endl;
  }
  fisse2.close();
*/
}

// ***********************************************************************
// This function checks if the file infilename has the format
//		h(int) k(int) l(int) Int(double) Sigma(double) Counter(long) ...
//		CorToNext(double) CorToSecond(double) etc. 	
// And returns the number of lines and the number of doubles in infilename through
// its first two arguments
// ***********************************************************************

void Diffractor::get_information_infile(long *number_of_lines_ptr, long *num_extra_lines_ptr, 
							long *highest_multiplet_ptr, char infilename[])
{
	long highest_multiplet_dummy;   // used in determining highest_multiplet

	FILE *infile_fp;		// where the file infilename is kept
	
  int dummy_h, dummy_k, dummy_l;		// used to read in numbers from the
	double dummy_int, dummy_sigma;		// file infilename
	long dummy_counter;					//	
		
	int previous_counter;				// used to check for doubles
	
	char line[5000];						// to read in first line of infilename
	
	
	// open infilename
	
	infile_fp = Fopen(infilename,"r");

	
	// read in first line	
	
	if (fgets(line, 5000, infile_fp) == NULL) 
	{
		printf("\nError reading first line of %s\n", infilename);

		exit(1);
	}



	// checks that the format of the input file is not completely crazy

	if ( 6 == sscanf(line, "%d%d%d%lf%lf%ld", &dummy_h, &dummy_k, &dummy_l, 
		&dummy_int, &dummy_sigma, &dummy_counter) )
	{

	}
	else
	{
		printf("The input format of %s is invalid\n", infilename);
		printf("Please see manual.doc explaining the required format of the input file\n");

		exit(1);
	}
	
	
	// move the file pointer associated with infile_fp back to the beginning of the file
	
	fseek(infile_fp, 0L, SEEK_SET);
	
	
	// initiate previous_counter to dummy number
	
	previous_counter = -1;
	
	
	// initiate number of lines and multiplet
	
	*number_of_lines_ptr = 0;
	*num_extra_lines_ptr = 0;
	*highest_multiplet_ptr = 0;
	highest_multiplet_dummy = 0;
	
	// find the number of lines etc
	
	while ( fgets(line, 5000, infile_fp) != NULL )
	{
		if (count_number_in_line(line) == 0)
		{
			break;
		}
	
		// increase number of lines by one
		
		++(*number_of_lines_ptr);
		
		
		// check for the presence of a multiplet
		
		sscanf(line, "%d%d%d%lf%lf%ld", &dummy_h, &dummy_k, &dummy_l, 
			&dummy_int, &dummy_sigma, &dummy_counter);
		
		if (previous_counter == dummy_counter)
		{
			// increase number of doublets by one
			
			++(highest_multiplet_dummy);
			++(*num_extra_lines_ptr);
		}
		else
		{
			if (highest_multiplet_dummy > *highest_multiplet_ptr)
				*highest_multiplet_ptr = highest_multiplet_dummy;

			// put highest_multiplet_dummy back to zero again

			highest_multiplet_dummy = 0;

			// set previous_counter to current counter
		
			previous_counter = dummy_counter;
		}
	}
    

	// close file
	
	fclose(infile_fp);
}


FILE* Diffractor::Fopen(char *fname, char *fmode) {
  FILE *fp;
  
  if ((fp=fopen(fname, fmode)) == NULL)
  {
    cout << "Couldn't open file " << fname << " in mode " << fmode << endl;
  }
  return(fp);
}


///////////////////////////////////
// Read in a 'new' type Dash hcv file.
//
// However split any multiplet into singlet. 
// Assume multiplet have the format:
//     h k l ..... 32 100 4
//     h k l ..... 32 4
///////////////////////////////////
void Diffractor::populate_store_data(char* DiffractorFile, vector<int>& num_cor_element)
{		
	int counter;			// the counter in the data file

	long num_extra_lines;				// number of extra lines from multiplets

	long highest_multiplet;    // highest multiplet in input data file
														 // highest_multiplet 0 = only singlets
														 // highest_multiplet 1 doublet and so on

	long number_of_lines;			// equal to num_refined_reflections + num_extra_lines

  long PAWLEY_COR_ELEMENTS = 20;


	// checking and getting number of refined reflections etc from infile
	
	get_information_infile(&number_of_lines, &num_extra_lines, &highest_multiplet, 
		DiffractorFile);
		

	// open infile which will be read below 
	
	ifstream data_file(DiffractorFile); 
	
	
	data.N = number_of_lines;   

  num_cor_element.resize(data.N+1); // +1 to conform to Numerical recipes standard

	newM(data.inv_cov, data.N, PAWLEY_COR_ELEMENTS);

  {
  for (int i = 1; i <= data.N; i++)
    for (int j = 1; j <= PAWLEY_COR_ELEMENTS; j++)
      data.inv_cov[i][j] = 0.0; // done since the code calling this method make use of none populated elements are zero
  } // somehow microsoft visual studio gets confused here and don't properly clear i and j

  newA(data.ave_intensity, data.N);
	
	newA(data.h, data.N);
	newA(data.k, data.N);
	newA(data.l, data.N);


	// initialize last_counter to dummy number


	int i = 0;   // pointing at first peak in infilename

  string aLine;

	while( getline(data_file, aLine, '\n') )
  {
    int number = count_number_in_line(aLine.c_str());

    if (number == 0)
      break;

    i++;

    num_cor_element[i] = number - 6; // h,k,l,I,sigma,count not counted


    stringstream ins(aLine);

		// read in h,k,l
		
		ins >> data.h[i] >> data.k[i] >> data.l[i];
		
		ins >> data.ave_intensity[i] >> data.inv_cov[i][1];


    // read in the data counter
		
		ins >> counter;
		
		
			
		for(int j=2;j<=num_cor_element[i]+1;j++)
    {
			ins >> data.inv_cov[i][j];

      // In principle the input file should not contain -100 or 100 correlation 
      // elements - since if two peaks are 100% correlated than better to treat
      // these as a one refineable. Further to avoid rule out the possible that
      // a matrix inversion struggle then do the stuff below

      if (data.inv_cov[i][j] < -100.0001 || data.inv_cov[i][j] > 100.0001)
      {
        printf("\n\n***********************************************************\n\n");
			  printf("ERROR IN INPUT FILE:\n");	
			  printf("  Correlation element for reflection %d %d %d\n", 
          data.h[i], data.k[i], data.l[i]);
			  printf("  is larger than 100 percent, which is clearly nonsense\n.");
			  printf("***********************************************************\n\n");

			  exit(1);
      }
    }
  }

	// double check that infilename has been stored correctly in data[0]
	
	if (i != data.N)
	{
		printf("Error in reading %s\n",
			DiffractorFile);
		exit(1);
	}

}


// return number of items in line
int Diffractor::count_number_in_line(const char* line)
{
		int i = 0; // plays the role of stepping through the string line
		int j = 0; // plays the role of counting the number of numbers in line
		do
		{
			i++;
			if ((line[i] == ' ' || line[i] == '\t') && (line[i-1] != ' ' && line[i-1] != '\t'))
			{
				j++;
			}
			if (line[i] == '\0') break;
		} while (line[i] != '\n');
		if (line[i-1] != ' ' && line[i-1] != '\t') j++; 
		
    return j;
}


